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ABSTRACT 

In gasdynamic systems, information travels in one direction 
for supersonic flow and in both directions for subsonic flow. A 
shock occurs at the transition from supersonic to subsonic flow. 
Thus, to simulate these systems, any simulation method implemented 
for the quasi-one-dimensional Euler equations must have the ability 
to capture the shock. In this paper, a technique combining both 
backward and central differencing is presented. The equations are 
subsequently linearized about an operating point and formulated 
into a linear state space model. After proper implementation of 
the boundary conditions, the model order is reduced from 123 to 
less than 10 using the Schur method of balancing. Simulation 
results comparing frequency and step responses of the reduced order 
models and the original system model are presented. This paper 
essentially follows the approach of Chicatelli, 1990, and 
Chicatelli and Hartley, 1990, while using an alternative for the 
flux splitting method. 


INTRODUCTION 

The digital simulation of gasdynamic systems is typically 
complex and computer intensive. Furthermore, gasdynamic systems 
have general characteristics which must be accounted for when 
choosing a simulation method. One major characteristic to be 
considered is that in gasdynamic systems, there are regions of 
supersonic flow and subsonic flow. In a one-dimensional problem, 
regions of supersonic flow will exhibit travelling of information 
in one direction. On the other hand, regions of subsonic flow will 
have information travelling in two opposing directions. The 
discontinuity occurring at this transition from supersonic flow to 
subsonic flow is called a shock. A shock generally appears in any 
physical system where there is a substantial transfer of energy 
from one form to another. In these systems, shock position is 
usually the desired control variable because of its physical 
significance. Consequently, any simulation method used must have 
the ability to accurately track the shock position. This specific 
consideration will be discussed later in the paper, but first 
consider a general discussion on the simulation of gasdynamic 
systems . 



A popular simulation approach for partial differential 
equations is the finite differencing technique. However, spatial 
differencing on distributed partial differential equations requires 
some general knowledge of system behavior. Moreover, the effects 
of the differencing method on the simulation must be considered. 
Since forward or backward differencing allows information to flow 
in only one direction, either method leads to numerical 
instabilities for systems containing subsonic flow. Thus, forward 
and backward differencing are, by themselves, not suitable 
simulation techniques for gasdynamic systems. On the other hand, 
although central differencing permits the flow of information in 
both directions it typically leads to unstable difference equations 
and/or simulations corrupted with high frequency spurious noise 
[3]. The foregoing suggests that, to simulate gasdynamic systems, 
one may consider a method which combines forward or backward 
differencing with central differencing. However, it should be born 
in mind that, assumptions about information flow direction must be 
made. 


The simulation method implemented in this paper is essentially 
a modification of a method developed by Courant, Isaacson, and Rees 
[4] and Roache [7]. This method considers the actual physics of 
the gasdynamic process when performing the spatial differencing. 
Mass flow and energy are assumed to propagate signals downstream. 
Therefore, only backward differencing is used in these two 
equations. Any term associated with system pressure is assumed to 
communicate information in both directions. Thus, pressure terms 
are estimated using central differencing. It turns out that this 
method has the ability to capture the shock, remain stable, and 
provide accurate results [3], Since it performs the differencing 
based on fundamental physics, the method is referred to as physical 
lumping. This technique is applied to the general quasi-one- 
dimensional gasdynamic equations for density, mass flow, and 
energy. Once this spatial differencing scheme is applied, the 
gasdynamic equations are subsequently linearized about a steady 
state operating point. 

The specific system under consideration is the NASA Lewis 40- 
60 Inlet [8]. This physical system may be represented by 41 
spatial lumps approximately 0.1427 feet apart. Considering that 
there are three governing equations for each spatial lump, the 
overall system is 123rd order. Since the dynamics at each lump are 
only a function of the previous lump and the next lump, the 
structure of the model lends itself to a tridiagonal state space 
formulation. This high order state space representation is then 
reduced using the Schur method of balancing [5]. 

The Schur method of balancing and its use in model reduction 
were first presented in [5]. In this truncation based model 
reduction method, the first concern is the size of the 
characteristic Hankel singular values (HSV) . Typically, any large 
break in the HSVs, usually taken as a 10 to 1 ratio, is a feasible 
position to truncate the model in balanced coordinates. The 
resulting reduced order model (ROM) must then be studied to ensure 



that all desired characteristics of the linear full order model 
(LFOM) and the nonlinear full order model (NLFOM) are retained. In 
this paper, 4th and 6th order linear ROMs are calculated. The step 
responses and frequency responses of these ROMs are considered. 


NASA Lewis 40-60 Inlet 

The starting point for the analysis begins with the three 
governing nonlinear quasi-one-dimensional gasdynamic equations for 
continuity, conservation of momentum, and conservation of energy. 
These equations, referred to as the Euler equations, are 

Continuity: 

d (PA) + d (puA) =M (1) 

at 8x 


Conservation of Momentum: 

a ( puA) + a [A (P + pu 2 ) ] = p 8A + F 

at 8x 8x 

and Conservation of Energy: 

8 (EA) , a[Au(E + P)] . _ D aA A A 

at ax at Vs 

with 



The major variables of concern are 

p = density 
u = velocity 
P s pressure 
E s energy 
A = area. 

The velocity may be expressed in terms of the mass flow and the 
density as 

U = (5) 

P 

Furthermore, the general gasdynamic equations may be written in 
terms of p, m, E, and P. The terms M s , F s , and Q s are the input 
source terms. For the specific problem considered here, inputs 
will be applied at the boundary conditions. After this 
simplification, the system equations may then be expressed by 
finite difference terms. 



The spatial differencing method implemented in this paper 
applies backward differencing on the mass flow and energy terms 
since they are considered to flow in only one direction. On the 
other hand, central differencing is applied to any pressure-related 
terms since they are considered to propagate in both directions. 
Thus, the equations for density, mass flow, and energy are first 
discretized in space as 
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The NASA Lewis 40-60 inlet may be represented by 41 spatial lumps 
with a spatial step, H - 0.1427 feet as in [1]. This in turn 
results in a 123rd order model. It should be noted that this 
specific method uses Euler's method to approximate the time 
derivatives. To expedite the analysis process, and to keep from 
adding a dynamic equation for pressure, the system equations must 
be further simplified. 

All of the pressure variables shown in (6) may be expressed in 
terms of the relationship shown in (4) . Subsequently, (6) may then 
be rewritten as follows 
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(7) 

Equation (7) can now be implemented as a shock capturing 
computational fluid dynamics (CFD) scheme for the quasi-one- 
dimensional Euler equations. The resulting model can be used for 
control evaluation in a real-time simulation and thus for control 
system design. 

The next step is to linearize equation (7) about a steady 
state operating point in order to develop a model for controller 
design. However, it is necessary to first consider the structural 
properties of (7) . Recall that the 40-60 inlet is comprised of 41 



spatial lumps with three governing equations for each lump, namely, 
density, mass flow, and energy. The structure of these equations 
enables the construction of a large linear state space model. 
Suppose a state vector is defined as 

X T = [Pi 1% E x : p 2 m 2 E 2 : . . . : P 41 ITl 41 ^*4l] ” (®) 

Following this definition of x, the system may be put into the 
general state space form 

x = Ax + Bu / 9 \ 

y = Cx. 


Since the dynamics of the system are a function of the present 
state and the states to the left and right, the A matrix has a 
tridiagonal form and may be constructed as 


■J x Q 2 0 0 . . . 0 

P 1 j 2 q 3 0 . . . 0 

0 P 2 J 3 Q 4 . . . 0 


0 0... p n _ 2 J n _ x Q n 

0 0 0.. . P n -i J n . 


( 10 ) 


Notice that the first and last rows do not include a P ; and 
term respectively. The effect of these matrices must be considered 
when applying any boundary conditions to the system as will be seen 
later in the paper. Taken directly as the small perturbation 
linearization of the discrete space equations given in (7), the P i , 
J. , and Qj matrices shown in (10) are 
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The measured output is the change in pressure, <5P, directly 
downstream of the shock. The position of the shock is directly 
related to the change in pressure at this position [6]. As the 
steady state shock position in the 40-60 inlet is located around 
the 24th lump, the C matrix for the output 6P may be constructed as 


C = 


0 0 


o ° ■ 2m » 0.4 0 
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P 24 


P 24 



( 14 ) 


The input to the system will be a change in pressure at the 
farthest point downstream corresponding to the last physical lump. 
This input is introduced into the system model by the B matrix. 
However, the effects of the boundary conditions should be 
considered first. Boundary conditions for this problem must take 
care of introducing a reflection in pressure information at the 
last lump. Consequently, the last row of the A matrix must be 
changed. The last row of the small perturbation model, with 
corresponding subscripts for elements in P and Q, is given as 

5E 4 i = ( ^40 ) 3 , 1^ P 40 + (^40^3,2^ m 40 + ^40^3,3^^40 (15) 

+ ^ 41 ^ 3 , 1 ^ P 4 I + 43 .) 3,2 6 m 4 i + ( J 41 ) 3 , 3 ^E 41 , 
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The last row in the (J 41 ) rou col submatrix, with designated subscripts, 
must be modified to create a new (J A1 ) rou<col submatrix, namely 
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( 17 ) 


Under further simplification, the equations in (17) become 
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These modifications must be made in the last row of the A matrix 
for proper implementation of the boundary conditions. 


Since <SP 41 is the input point, the $P 41 term of the <SE 41 
equation in (11) may be absorbed into the B matrix, viz. 
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The nonzero term in (19) is simply the coefficient on 6P given in 
(16) and multiplied by (J 41 ) 33 . The state space model for the 
system is now completely defined. In the next section, the state 
space model is balanced, using the Schur method of balancing, and 
truncated. The resulting ROMs are then studied using the step and 
frequency responses. 


Model Reduction 

Obtaining an accurate reduced order model is very important 
for controller design, particularly for reducing the complexity of 
the resulting controller. No reliable techniques are currently 
available for reducing the order of the original nonlinear system 
while preserving large perturbation information. However, many 
methods are available for reducing the order of linear systems. 
Among these methods, the Schur method is a robust and well 
conditioned method to reduce the large state space models of 
gasdynamic systems. A reduced order model may be found directly by 
using projections defined by the left and right eigenspaces of the 
large eigenvalues of the product of the observability and 
controllability Grammians. Since the full order models are 
typically numerically ill-conditioned, some type of scaling is 
performed prior to model order reduction. A first consideration in 
choosing the ROM order is to view the largest HSVs of the system. 



An abbreviated table of the 13 largest HSVs of the linear 
model of the NASA Lewis 40-60 inlet is shown in Figure 1. 


order 

HSV 

ratio 

1 

1. 0603e+01 


2 

4 . 6175e+00 

2.2962 

3 

1. 9535e+00 

2.3637 

4 

1 . 0341e+00 

1.8891 

5 

4 . 9687e-01 

2.0812 

6 

2 . 6303e-01 

1.8890 

7 

9 . 2066e-02 

2.8570 

8 

9. 1363e-02 

1.0077 

9 

2 . 1703e-02 

4.2097 

10 

1. 1679e-02 

1.2926 

11 

1 . 6844e-03 

6.9336 

12 

1 . 5313e-03 

1.1000 

13 

1 . 5189e-04 

10.0816 


Figure 1. - Largest Hankel Singular Values of the NASA 
Lewis 40-60 Linear Inlet Model. 


The ratio column in Figure 1 is simply the ratio of the HSV to the 
left of the number divided by the HSV above that number. A ROM is 
typically found by truncating the LFOM where there is a 10 to 1 
break in the HSVs. At first glance the best break appears to be 
for a 12th order model. Other possible considerations include 8th 
and 10th order ROMs. It turns out that all of these ROMs trace the 
step response so close that they are not discernable on the graph. 
Subsequently, two other ROMs were calculated of order 4 and 6. 
Since the output is known to have a delay in its time response, the 
order of the model cannot be reduced much more than this. Even the 
4th order model exhibits some oscillatory behavior when trying to 
represent the time delay. The step response for these two ROMs is 
shown in Figure 2. Also included in figure 2 are the step 
responses for the LFOM and the NLFOM . The 6th order model traces 
the actual LFOM response very close. As mentioned previously, the 
4th order response is somewhat oscillatory during the time delay 
but still gives a reasonable approximation to the actual step 
response. For completeness, the frequency response is shown in 
Figure 3 . 
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Figure 2 - Step Response of LFOM, NLFOM, and ROMs 
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Figure 3 - Frequency Response for LFOM and ROMs 




From Figure 3, it can be seen that the original linear system is 
very low pass. Notice that the frequency responses separate only 
when the system starts to attenuate. If higher order ROMs are 
considered, the only effect on the frequency response is that more 
roll-off is preserved by keeping more poles. Overall, the accuracy 
of the ROMs appear to be acceptable considering the amount of order 
reduction attempted. R0M4 is shown in Figure 4 while R0M6 is shown 
in Figure 5. 
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Figure 4 . - 4th order reduced order model for the 
NASA Lewis 40-60 Inlet 
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1 . 29e+2s 5 - 9 . 44e+5s 4 + 3.59e+9s 3 - 8.12e+12s 2 

+ 1.06e+16s - 6 . 36e+18 


s 6 + 1 . 90e+3s 5 + 6.54e+6s 4 + 7.22e+9s 3 + 7.77e+12s 2 

+ 3 . 12e+15s + 4 . 45e+17 

Figure 5. - 6th order reduced order model for the 
NASA Lewis 40-60 Inlet 



Conclusions 


The given physical lumping method of differencing proved to be 
feasible for representing the dynamics of the NASA Lewis 40-60 
inlet. An advantage of this differencing approach is that it 
readily allows the study of nonlinear model reduction methods, as 
the states are immediately available. This is not true of most 
other CFD methods. Furthermore, physical lumping is more 
intuitive. It is essentially a straight forward differencing 
approach which requires less up-front calculation than the split 
flux method presented in [1], The resulting ROMs not only turned 
out to be of smaller order but they also efficiently captured the 
dynamics of the system. The Schur method of balancing proved to be 
a good choice for a model reduction scheme resulting in a 
substantial reduction in order from 123 to order 6 or smaller. 
Depending on the intended use of the reduced order models, the 6th 
order model appears practical for most purposes. The 4th order 
model may be used if time delay information is not important. 
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